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COMPUTER SIMULATION OF PLASMA AND N-BODY PROBLEMS 


By 

Wynford L. Harries^, John B. Miller^, and Christopher Costner^ 

INTRODUCTION 

In recent years, large N-hody computer simulations (Miller and Prendergast, 
1968; Hohl and Hockney, 1969) have become an important tool in investigating 
the structure of spiral galaxies, especially in determining the development of 
large-scale instabilities resulting in spiral and bar formation. Until recently, 
most of these simulations used essentially two-dimensional models with the 
"stars" confined to the plane of the galactic disk (Miller, Prendergast, and 
Quirk, 1970; Hohl, 1971). These simulations have sho™ that the disks of stars 
have a tendency for the development of fast growing nonaxisymmetric instabilities 
resulting in bar formation. The bar instabilities occur even for velocity 
dispersions that are considerably larger than those found in the solar neighbor- 
hood or those predicted by Toomre (1964) as being locally stabilizing. Because 
of the difficulties in solving the highly nonlinear problem, global instability 
studies of disks of stars have been primarily numerical. Some limited work 
has been done for uniformly rotating disks (Hunter, 1963; Kalnajs, 1972), but 
generally linear stability analyses were used in the studies of disks of stars. 

Any spiral structure in computer-generated galaxies is generally short 
lived and the final state is a rotating bar. The bar thus obtained rotates 
more slowly than the stars. For one case investigated by Hohl (1971), the 
bar rotates at 2.25 t and the stars rotate at 1.5 t where x is the rotational 
period of the initial disk. 
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^ Graduate Research Assistant, Department of Physics and Geophysical Sciences, 
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It has been argued that core/halo components have a stabilizing 
effect on galaxies and result in longer lived spiral structure fOstriker 
and Peebles, 1973). However, numerical experiments with large fixed 
stellar components representing the core/halo component (Hohl, 1970; 

Hockney and Brownrigg, 1974) show that multiarmed spiral structure devel- 
ops and persists for many rotations but only in an evolving manner- That 
is, the spiral structure is either wound up into a tight pattern or it is 
wound up and then reappears again. A recent study of the effect of fixed 
core/halo components (Hohl, 1976) does show that the bar instability is 
indeed inhibited by a sufficiently large fixed component. 

The purpose of the present study is twofold.' First, we want to deter- 
mine the effect of a self-consistent (rather than fixed) core/halo compon- 
ent. This will show whether there are any instabilities (such as ’’two- 
stream") or other important interactions present that may be suppressed 
with a fixed core. Second, we want to determine the effects of finite 
thickness of ‘the disk and of three-dimensional essentially spherical 
core/halo components. 


MODEL 

The model used for the present galaxy simulations consists of 100,000 
representative stars that move inside an array of cells. For the disk 
simulations the stars are confined to move in the plane of the disk repre- 
sented by a 64 X 64 active array. In the three-dimensional simulation the 
stars move inside a 64 x 54 x 16 array of cells. The sum of the stars in- 
side each cell defines the mass density at the center of each cell. Fast 
Fourier transform methods are used to obtain the gravitational field at 
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the center of each cell for a given density distribution. The force act- 
ing on a particular star is determined by bilinear (or trilinear) 
interpolation from the values of the gravitational fields at the surround- 
ing 4 (or 8} cell centers. After the force acting on a star is determined 
it is advanced by a small timestep, the new density is recalculated and 
the process is continued until the desired evolution is achieved. If a 
star leaves the array of cells, approximate methods are used to determine 
the force acting on the star. Details of the disk model are described in 
detail by Hohl and Hockney (1969:i and by Hohl (1970:i , The extension of 
the model to three dimensions is described in the appendix. 

RESULTS 

Observational evidence (deVoucouleurs 1959; Freeman 1970; Korraendy 
1977) indicates that the luminosity (and presumably the density) in the 
outer regions of many spiral and SO galaxies decreases exponentially with 
radius. Also, previous simulations (Hohl 1971) showed that intially un- 
stable stellar disks evolved into stable systems with radial density 
variations that closely approximated the sum of two exponentials. The 
inner exponential with a scale ‘length of about 1 kpc describes the non- 
or slowly-rotating spheroidal or core component and the remaining 
exponential with a scale length of about 8 kpc describes the extended disk 
population. Thus, it seems reasonable to use an exponential density 
variation for the disk of the present computer simulations. Similarly, 
the central core used is described by an exponential density variation. 

Figure 1 illustrates the evolution of a disk of 100,000 stars with an 
initially exponential surface density distribution y(r) = with 

a cutoff at r = 10 kpc. The initial angular velocity of the disk was 
obtained from 



J = t, 2 ^ __l 1 _ a ^(r)] + [a - <^Q^(r)] ( 1 ) 

rjiCr) 3r ^ r ® 

with 

~ g M ■ (2) 

2o 3 (r) 

0 

Here, “qW is the angular velocity required to balance the cold (zero 
velocity dispersion) ui(t) is the actual angular velocity, and 

K(r) is the epicyclic frequency. The initial value of the radial veloci- 
ty dispersion was taken to be that determined by Toomre (1964) as 

the minimum required to stabilize all axisymmetric instabilities, 

c^ 3 .Cr) = - 3.36 GpCr)/K(r) (3) 

of the cold disk at a 
radius of 5 kpc, that is half way to the edge of the initial disk. 

As expected (Hohl 1970, 1971), only the small-scale instabilities are 
prevented by a = a . and the system quickly forms a two-arm spiral 
which eventually tends to evolve into a rotating bar. The evolution of the' 
azimuthally averaged radial density variation for this system is shoim in 
Fig. 2. As previously observed (Hohl 1970, 1971) the eventual density 
variation approaches one which can be closely approximated by the sxam of 
two exponentials. One exponential describing the central core component 
and the other describing the extended disk. The evolution of the radial 
velocity is such that there is some heating near the center, and a consider- 
able increase in the velocity dispersion for stars expanding into the 
extended disk component. Numerous other diagnostics have been performed 
on the system. For example. Fig.- 3 shows the time evolution of the moment 
of inertia I divided by the moment of inertia at t = 0, and a similar 
ratio for the angular momentum. P. As can be seen, P is conserved in the 


The time t is given in rotational periods 
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simulation but I is still increasing at a near linear rate after three 

rotations. The evolution of various components of the total kinetic 

energy divided by the total potential energy is shown in Fig. 4. The 

components T^ and Tg represent the kinetic energies due to the 

2 '> 

velocity dispersions Z m.a. and S m.o. respectively, while T . 

i ^ ^0 • i 

is the kinetic energy of rotation. Note that the ratio of the kinetic 
energy in rotation to the absolute value of the total gravitational energy 
of the system is approaching the value 0.14 predicted by Ostriker and 
Peebles (1973] for stability. At the same time, there occurs considerable 
heating of the system. 

One of the aims of the present study is to determine the effect of 

adding the third degree of freedom by allowing a finite thickness of the 

exponential disk. Using again an exponential projected surface density 
* /2 

variation p = e the stars are now distributed in the z-direction 

2 

according to one-dimensional distribution sech z/c where c is a param- 
eter determined from y ^ (Hohl 1967) . The central thickness of the disk 
is 2 kpc and the density is cut off at given by 



where R = 10 kpc is the radius of the disk. The radial and azimuthal 
velocity components are determined in a manner similar to that for the 
infinitesimally thin disk and the z-component of the velocity dispersion 
is determined by a force balance in the z-direction. Note also that all 
initial velocities are truncated such that stars have kinetic energies no 
greater than that which would allow them to reach the boundary of the 
system in the gravitational potential at t = 0. 
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Figure S shows a side view of the initial disk and the evolution for 
up to 3 rotations. Note the rapid expansion in the plane of the disk. 

result of the bar instability as shown in Fig. 6 which gives 
the evolution of the disk projected in the x-y plane. Note that the 
evolution is very similar to that shown in Fig. 1 for the infinitely thin 
disk. Similarly the evolution of the surface density variation and the 
increases in the moment of inertia are nearly identical to those shown in 
Figs. 2 and 3 for the thin disk. The ratio of the various kinetic energy 
components for the total potential energy are shown in Fig. 7 for the 
finite thickness disk. Note that again the evolution is similar to that 
for the infinitely thin disk as shown in Fig, 4. An additional variable, 
the z-component of the kinetic energy, is given in Fig. 7 and shows that 
since this component remains small compared to the others one would, expect 
little difference in the evolution of the finite thickness disk when com- 
pared to the infinitely thin disk. 

As shown in Figs. 1 and 6, exponential disks ^^ith velocity dispersion 
CQ ~ 1) are violently unstable to the bar- forming instability. Previous 
work CHohl, 1976) with a superimposed fixed (nonself-consistent) central 
mass distribution indicated a stabilizing effect toward the bar- forming 
instability. A more realistic simulation is to allow core-disk interaction, 
thus, presently we are interested in the stabilizing effects of a completely 
self-consistent core or "spheroid" component. Again, the effect is investi- 
gated for both the infinitesimally thin disk (two-dimensional) and for the 
three-dimensional disk. 
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For. the core-disk system, 50 percent of the mass (50,000 stars^ is 
contained in the nonrotating core and the remaining mass (50,000 stars) 
is contained in the disk. The disk component is again given the surface 
density variation whereas, the initial nonrotating core 

component is given a density variation 

disk and core density are cut off at r = 10 kpc and r = 3.5 kpc, 
respectively . The initial velocity dispersion and rotation of the disk 
is obtained by again using Eq-. (1), (2), and (3) with y = Similarly, 

as before, the z-dimensions of the disk are determined from Eq. (4) . The 
initial velocity dispersion of the nonrotating core was obtained by taking 
^ simply balancing the core in the presence of the disk. In 

order to assure that the core component was in a stable state at the start 
of the core-disk simulation, the core was allowed to evolve for several 
rotational periods (2ir/ti)^ at 5 kpc) with the disk component held fixed. 
Starting from these initial conditions, the system evolved as shown in 
Fig. 8. Note that even though a two-arm spiral structure still forms, the 
system as a whole evolves in a much less violet manner than that displayed 
in Fig. 1, This can also be seen in Fig. 9 which shows the evolution of 
the surface mass density for both the core and disk component. Note that 
with the exception of a slow outward diffusion of stars near the edge, the 
core remains essentially stationary, while the disk component displays the 
outward shift of mass generally associated with bar formation. Similar 
information is contained in Fig. 10 which displays the evolution of the 
radial velocity dispersion for the core and disk component. Note the sharp 
increase in the velocity dispersion at r = 2 kpc which is associated with 
a marked reduction in the angular momentum of the disk in this region,- 
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In general, the simulations show that the foarmation o£ bars or two-armed 
spirals results in moving angular momentum outward to larger radii. 

Fig. 11 shows a marked reduction in the rate of increase of the moment 
of inertia when compared to the disks without a central core component . 

The final system investigated is that of a three-dimensional 
exponential disk with a three-dimensional core or spheroid component. 

The spatial distribution of the stars for the disk component is obtained, 
as ll/as done for the disk shown in Fig. 5 and 6, except that now the disk 
contains only 50,000 stars. For the nonrotating central core the density 
is given by P = where ? = + y^ + with c = 5/7. 

The density is cut off at 5=7. Thus, the central core or spheroid has 
an axis ratio of 7:5. Again the Gaussian velocity dispersion for the core 
is. obtained by a simple balance of the self-gravity of the total system. 

The velocity .dispersion for the disk component is generated, as was done for 
the system shoim in Fig. 6. Before initiating the simulation of the com- 
bined core-disk system, the core was allowed to evolve for several rotations 
(free-fall periods) to assure that no instabilities or other problems 
associated with the core component were present. 

Figure 12 shows the evolution of the system perpendicular to the 
equatorial plane. Note the remarkable stability of the system when com- 
pared to the disk i^ithout the central core in Fig. 5. The evolution of 
the system in the equatorial plane is sho™ in Fig. 13 and displays the 
development of a comparatively weak two-arm spiral structure. It should 
be noted that because of the allowed initial relaxation, the core components 
of the two core- disk systems investigated here are expected to closely 
satisfy the collisionless Boltzmann equation. The same is not necessarily 
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true for the disk component since satisfying equation (1) only assures 

a balance of forces at t = 0. Also, we know that for a stellar disk 

a = a . does not assure stabilization of global nonaxisymraetric 
r r,min 

instabilities (Toomre, 1974; 1977). However, since one would hardly 
expect nature to generate a galaxy initially in an exact stable station- 
ary state, and since we are interested in the further development of 
instabilities and the final state toward which the system evolves, an 
exact stationary and stable initial state is not necessary. 

The evolution of the azimuthally averaged projected surface mass 
density for the three-dimensional core-disk system is shown in Fig. 14 
and is nearly identical to that of the two-dimensional core- disk system 
shown in Fig. 9. Note that there is very little change in the density 

for the core with the exception of a slight outward diffusion near the 

» 

edge. Azimuthally averaged values of the total density variation 
in the z-direction are shown in Fig. 15 for various values of r. Some 
of the fluctuations shown may be due to the relatively small sampling .volume 
used. If we look at the evolution of the radial velocity dispersion shown 
in Fig. 16 we see that (as expected] the velocity dispersion for the 
two-dimensional core C^ig. 10) is higher. Also, the large increase in the 
velocity dispersion of the disk near r = 2 kpc does not occur for the 
three-dimensional disk. Associated with this is the fact that there is 
very little change in the radial angular moment distribution during the 
evolution of the 3-D core-disk system, Avhereas, considerable outward shift 
of angular momentum occurs for the 2-D core-disk system. These results 
indicate that the global bar instability is much weaker for the 3-D 
system as for the 2-D sys.tem, as can be seen by comparing Figs. 8 and 13. 
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The time evolution of the various kinetic energy ratios for the 3-D 
disk-core system is sho™- in Fig. 17. As can he seen, there is little . 
change in the value of the various components during the evolution. Note 
that the value of the ratio of the kinetic energy in rotation to the total 
potential energy of the system is slightly higher than the value of the 
0.14 predicted for stability by Ostriker and Peebles (197*3). Also, the 
moment of inertia increases by only about one third of that shown in 
Fig. 11 for the 2-D system. As was the case for all four systems investi- 
gated, the angular momentm was conserved to a sufficient degree of 
accuracy. 
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Figure 1.- Evolution of an initially balanced, infinitesimally thin disk 
of 100,000 stars with an exponential radial density variation. The 
stars have an initial density variation given by Toomre’s criterion. 
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Figure 3.- Time evolution of the angular momentum (P] and the moment of inertia (l) 
for the unstable disk shown in figure 1. Note the rapid increase in I as the 
bar begins to form at t a i. 
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Figure 4.- Time evolution of various kinetic to total potential energy ratios. Note 
that the ratio of rotational to iDotential energy is approaching the value of 0.14 
predicted by Ostriker and Peebles as required for stability. 
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Figure 5.- Side view showing the evolution of the three-dimensional exponential disk. 
The bar instability results in a rapid expansion in the plane of the disk. 
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Figure 6.- Evolution of an initially balanced three-dimensional stellar system 
of 100,000 stars with an exponential radial density variation. 
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Figure 7.- Time evolution of various kinetic to total potential energy 
ratios. Note the evolutions of the energy ratios are similar to 
those shown in figure 4. 
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Figure 8.- Evolution of an infinitesimally thin exponential disk 
with a self-consistent exponential core component. Note that 
the evolution if considerably less violent than that displayed 
in figure 1. 
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Figure 10.- Evolution of the azimuthally averaged radial velocity dispersion for 
the two-dinionsional exponential disk plus core system. 
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Figure 11.- Time variation of the moment of inertia and the angular 
two-dimensional exponential disk plus core system* 
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Figure 12.- Side view of the evolution of the three-dimensional exponential disk plus 
core system. Note the remarkable stability when compared to the three-dimensional 
disk-only system sliown in figure 5. 
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Figure 13.- Evolution of the three-dimensional disk plus core s/stem viewed 
in the equatorial (x-y) plane. Note the development of the comparably 
weak spiral structure. 
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APPENDIX 


COMPUTER PROGRAM EOR GENERATING THE THREE-DIMENSIONAL 
GRAVITATIONAL POTENTIAL DISTRIBUTION OF ISOLATED GALAXIES 


MATHEMATICAL SUMMARY 

The scaled gravitational potential at the center of ceil (x,y,z) is 
defined by the triple siuiunation over the three-dimensional array of cells 

(b — / ■ / , / • I y CAi^ 

1^0 j^O 1 -x,j-y,k-z 


,.2 . 2 , 2 ^ 

H. . , = (i +j +k ^ 

^0,0,0 " 


for i -5- j + k r 0, 


and p. ... is the mass density in cell (i,j,k). Because direct summation 

1,1 ,-k 

is much too time consuming to be practical, the triple summation is eval- 
uated by the convolution method using fast Fourier transforms (ref. (Al)) . 
That is, the Fourier transform of the potential equals the product of the 
Fourier transforms of p and H 

= A. . • (“5 

The gravitational potential (j> is obtained by taking the inverse 

x,y,z 

Fourier transform of equation (A2) . Rather than the usual complex Fourier 

series, here a real expansion is used. For example, the Fourier transform 

of the density P is given by 

Xj y j z 


C,n,c 


2h-l 2n-l 2n^l 

/ 7 7 c(x,n]c(y,n}c(z,h) p ,f (5,x,n)f (n,y,n)f (^,z,h) 

ht) ^ 

(A3) 
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where 


fCC,x,n3 =■ 


:os C? x/n3 , 0 ^ ^ n 

;in ['!r(C-n)x/n] , n < C 2n 


c(x,n) = II'/t. if X = 0 or x = n, 
c.Cx,n) = 1, 


othenvise, the symbols n and h define the n x n x h active array 

and also the (2n) x C2n.) x (2h3 larger array over which the Fourier 

transform must.be taken so that the potential for an isolated galaxy is 

obtained fsee fig. Al). Note that the density may be nonzero only in the 

smaller n x n x h array. Because of the s>Tunetry of H , the 

x,y j z 

'Vi 

Fourier transform can be obtained bv a finite cosine transform 




2 2 2 
>. c^(x,n) c^(y,n) c^(z,h) H. 


z=o y=u x= 


x,y,z 


cos (TTCx/n3 cos (imy/n) cos C^-?z/h3, 

0 <, x-n 

0 ^ ? ^h 


= ft = ft = ft 

C+n/1,5 ? + n,ri+n,^ 5+n,n,C h C+n,Ti+n,C+h 


= ft =ft =ft =F 

5,n+n,c+h C,n,?^h 

The next step in obtaining the potential xs to multiply P,, „ „ ‘by 
to obtain 

K,r\,Z 

r\j /Xj 'll 
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Computer Program Subroutine ’?hich Uses Only Core Storage 
Table A1 gives a Fortran listing of a con 5 )Uter program which may be 
used to obtain the potential by use of a [2n) x (2n) x h array of cell. 

N 

The variables 12A and 13A define the x,y and z dimensions, respectively, 
of the array used for the potential calculations, ^^hen the subroutine 
GETPHI is called, RHOCI,J,K) contains the mass density and GETPHI places 
the values of the corresponding gravitational potential in RHOCI,J,K]. 

The subroutine FTRANS(I, I2B]) has been written by R. Hockney (ref. A2) 
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and it performs a finite Fourier analysis or synthesis on the common 
input array and places the result in the common output array Y. 

The subroutine performs a cosine analysis for I = 2-, a periodic analysis 
for 1=3, and a periodic synthesis for 1=4. The subroutine 
GETSET(I,I2B^ initializes FTRANS and is called every time the arguments 
of FTRANS (I, I2B) are changed. The Fourier transform H i-s calcu- 

lated on an x (n+1) x (h+1) array only the first time that the 

subroutine is called and is kept in storage for subsequent use. 


The Fourier transform of p in the x-direction is generated 

x,y,z 

by obtaining the partial transform p_ for 0 <. 2n-l, 

%>y 

i\j 

0 ^ y ^ n-1 and 0 ^ z ^ h-1. p_ is zero outside of 

§ > y > 2 

this region because p is nonzero only over the 

X ,y , z. 

n X n X h active array. Next, the Fourier transform of 

p is performed in the y-direction obtaining the x-y partial trans- 

form p„ , for 0 ^ C ^ 2n-l, 0 <. n ^ 2n-l and 0 :< z < h-1. Since 
C,h,z 

p.. is zero for h < z < 2h-l, by use of one-dimensional arrays Y and 

- - 

Z the Fourier transform of p„ can be taken in the z-direction to 

S,ri,z 

obtain the total transform p.. „ for 0 <. ? <L 2h-l. Next, zs 

multiplied by 5^.. to obtain ^ and the inverse Fourier transform 

is performed in the z-direction. The resulting partial x-y transform 
% is placed in the 2n x 2n x h RHO(I,J,K^ array for 

0 C <:. 2n-l, 0 u .2. 2n-l and 0 ^ z h-1. with values for h <. z ^ 2h-l 
discarded. (The use of these one- dimensional arrays was first presented in 
reference A3 for a two-dimensional potential solver). Next, the inverse 
Fourier transform of 4) is generated in the y-direction by obtaining 

9 j R > z 
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for 0 ^ C 2n-l, 0 ± y ^ n-1 and 


•\i 

the x-partial transform ^ 

0 i z i h-1. The final step is to perform the inverse Fourier transform 

in the x-direction for 0 :< y £ n-1 and 0 :< z ^ h-1 to yield, the 

correct gravitational potential <fi for an isolated galaxy over 

x,y , z 

the n x n X h array. 

Overlayed Computer Program WTiich Uses 
Core and Disk Storage 

The use of the listing of Table A1 with the 64 x 64 x 16 active 
density/potential array used in this paper would have necessitated the 
dimensioning of tlie RHO array at 128 x 128 x 16 and the H array at 
65 X 65 X 17. As such, large dimensions would have excluded use of the 
CDC 6600 computer, the listing of Table A1 was modified to include use of 
overlayed programs and disk storage resulting in a maximum core storage 
at any one time of array elements equaling about five fourths of the 
active array. The listing of this program in Table A2 includes (a) a 
section of an initializing overlay in which relevant constants are computed 
(b) a section of the star advancing overlay in which "chunks" of the den- 
sity array are written on appropriate disk files, another section of 
the star advancing overlay in which "chunks" of the computed potential 
array are read from disk files, (d) the GETH overlay which computer 
and (e) the GETPHI overlay which computes the potential array from the 
density array. 

The method used is the alignment in the direction of transformation of 
four identical arrays named "RHOl, RH02, RH03, and RH04, each of which is 
dimensioned Cn./2) x Cn/2) x h within the GETPHI overlay. (See figs. A2 
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and A5, For clarity, figures A1 through A6 are dra^vn for an active array 
dimensioned nxnxh=8x8x'4; table A3 compares the array dimensions 
of these figures and the listing of table A2.) The active array is 
dimensioned as the PHI array within the initializing and star advancing 
overlays (see figures A1 and A3) but is not dimensioned within the 
GETPHI overlay. As figure A2 suggests, the "chunks" RHOl, RH02, RH03 
and RH04 may be visualized as forming either a row or a column of the 
lower half (0 <. z ^ h-1) - of the extended array. Switching the lineup tO' 
a different row or column is accomplished by storing the array associated 
with each "chunk" location on a separate file; these eight files are also 
indicated in figure A2, 

As shora in figure A3 one "chunk” size array named 01 is dimensioned 
in the initializing and star advancing overlays. "Chunks" of the active 
array are transferred between the PHI array of these overlays and the 
arrays RHOl, RH02, RH03 and RH04 of the GETPHI overlay via "do loop" trans- 
fer to/from the _ 0K_ array and storage on files 1, 2, 5 and 6. 

At the beginning of a program run, the GETH overlay computes H in the 

(n+l) X (n+1) X (h+1) H array in the same manner as the listing of table 

Al. All of H , except for two boundary planes of elements (? = n, 
h > ? 

O^C^h and 0 ^ C ^ n, n = n, 0 <. C ^ h) , is then transferred 
in portions via "do loop" to the (n/2) x (n/2) x (h+1) HH array from 
which it is written on disk file 9 (see figure A4) . Elements of one 
boundary plane of H (5 = n, 0 .< n ^ n, 0 < ^ < h) are transferred to 

the (n+1) X (h+1) HN21 array which is in common with the GETPHI overlay; 
the transpose of that boundary plane is equal to the other boundary 
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plane ri = n, 0^C<.h) due to the symmetry of H across 

the ?=n diagonal plane.' During- each potential solution the portions of 
H on file 9 are read sequentially into an i‘n/2') x (n/23 x Ch+1) HH 

r\j 

array of the GETPHI overlay from which H elements, along with those in 
the HN21 array, are multiplied with p . This sequence (listed in 
table A4) utilizes the symmetry and periodicity of H (equation (A4)) 

a, 

to provide a full set of (2n) x (2n] x (2h) H elements to the GETPHI 
overlay in a manner which minimizes the reading of file 9. 

The GETPHI overlay consists of subroutines ANLX(JCOLUMN) , ANLSYN(IROW) 
and SYNX(JCOLUMN) which dimension in common the arrays HH, HN21, RHOl, 


RH02, RH03 and RH04 as pictures in figure 5. Figure 6 indicates the 
lineup of ”chunhs’' associated with each call to a subroutine. The potential 
solution is mathematically identical with that described for the listing 
of table Al. Calling ANLX(l) and ANLX(2) performs the Fourier transform 
of p in the x-direction to form p . Calling AMLSYN(13, 

y fZ Z yY yZ 

ANLSYN(3), ANLSYN(2) and ANLSYN(4) in sequence performs the following: 


(a] a Fourier transform of p in the y- and z-directions to form 

A* ^ . 

p^ ; (b) multiplication with H to form (p ; and (c) the in- 

verse Fourier transform of 6,. in the z- and y-directions to form 

A/ , 

(J) . Calling SYNX(l) and SYNX(2) performs the inverse Fourier trans- 

form of ^ ^ in the x-direction to form GETPHI overlay 

§,y,z x,y,z 


is outlined in more detail in table A5. 
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Efficiencies of the Two Computer Programs 
The program of table A2 is considerably more efficient than that of 
able A1 because the addition of some peripheral processing time and a 
nail increase in central processing time is much more than compensated 
or by a 75 percent decrease in the required core storage. The maximum 
umber of active array elements dimensionable on the CDC 6600 with the 
rograms of table A1 and A2 are respectively 16384 (e.g. 32 x 32 x 16) 
id 65536 (e.g. 64 x 64 x 16); the latter program can have other 
otentially useful active array dimensions of 32 x 32 x 8, 32 x 32 x 16, 
ad 32 X 32 X 32. Solution of the 64 x 64 x 16 active array by the 
DC 6600 requires about 300 (octal) words of core storage and with H 
Iready computed takes about 75 seconds of central processing time. 
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SUBROUTINE FOR CALCULATING THE THREE-DIMENSIONAL 
GRAVITATIONAL POTENTIAL USING ONLY CORE STORAGE 



subroutine GETPHJ 

COMMON Z( 1 025) « Y( 1025) «RHO ( 64 1 6 ) • 12A . I3A • I TEST 

DIMENSION H(33«33U7) 

IF ( ITEST.EQ.O ) GO TO It 
ITEST=0 
128= I2A-1 
N=2»*I2A 
N02=N/2 
N21=N02+1 
133= I3A-I 
NH=2**I3A 
NH02=NH/2 
NH21 =NH02+1 
■ RNI=1 ./(N»N*NH) 

DO 1 K=1 .NH21 

DO 1 U=I *N2l 

DO 1 I=I»N21 

RI = (K-1 )* (K-1 )+( J-1 )♦( J-1 > + M-n*(I-l) 

IFIRI .LT.l . ) Rl = l . 

H£ I , J»K)=HNI/SQRT (R I ) 
a CONTINUE 

CALL GETSETC2« 12B ) 

DO 2 <=1 .NH21 
DO 2 J=I.N21 
DO 3 I = I « N2 1 

3 Z ( I )=H( 1 , J.K) 

CALL FTRAN3(2« I2B > 

DO 4 I = t tNEl 

4 H{ I , J.K )=Y( I > 

2 CONTINUE 

DO 5 K=! .NH21 
DO 5 I = I »N2l 
DO 6 J=I»N21 

6 2( J )=HI I .J.K) 

CALL FTHANSC2* I2B) 

DO 7 J=1 .N21 

7 H( I , J.K)=Yl J ) 

5 CONTINUE 

call GETSET(2. I3B) 

DO I 0 J = 1 «N21 
DO 10 1=1 »N21 
DO 8 K=1 .NH21 
S Z(K)=H( I .J»K) 

call FTRANS(2. I3B) 

DO 9 K=1 «NH21 
9 HI I .U»0=YCK) 

10 CONTINUE 

11 CONTINUE 
WRITE(6*43) 

43 FORMATdOH H(1*J.K)) 

DO 42 K=I *NH2l 
DO 42 J=1 tN21 . 

WR1TEI6*41 ) J»K 

WR!TEI6*40) (HI I • J»K ) • I = 1 »N21 )' 

41 F0RMATI14H I=l4N21 J=I3«SH K=I3) 

40 FORMAT (2H BE1S.8) 

42 CONT I NUE 

CALL GETSETI3.I2A) 

DO 14 K=l ,NH02 
DO 14 J=1 *N02 
DO 12 1=1 »N 

12 2(1 ) = RHO ( I 4 J tK ) 

CALL FTRANS(34 I2A) 

DO 13 I=l4N 

13 RHO( I 4 J4K )=Y ( I ) 

14 continue 


original ^ AGET is 

OF POOR QUALITY 



DO 17 K=1 1 NHO 2 
DO 17 1=1, N 
DO 15 J=1 ,N 

15 2{ J)=RHOt I , ) 

CALL FTHANS (3, I2A ) 

DO 16 J=1 ,N 

16 RHOU , J.K )=Y ( J ) 

17 CONTINUE 

DO 20 I =1 ,N 
DO 20 J=l ,N 
00 13 K=1 ,NH02 
Z(l<)=RHO( I , J,>< ) 
le Z(K-+NHO 3 )= 0 , 

CALL GETSET (3. I3A ) 
call FTRANSO, I3A ) 

IF ( 1 .GT.N2I .AND.J.LE.N2I 1 GO TO 22 
IF ( t .LE.N21 ,AN0.J.GT.N2! ) GO TO 2a 
IF ( I .GT,N21 ,aNd.u,gt.N 21 ) GO TO 26 
DO 19 K=l ,NH02 

Z(k:)=Y(K)*h( j ,j,k) 

19 2<K+NH02 ) =Y I<+nho2 >*H< I , J,K } 

Z(1 )=Y( 1 )*Ht 1 ,J,1 ) 

2(Nh 21 )=Y(NH21 )»H(r.J.NH2l ) 

GO TO 21 

22 DO 23 K=2,NH02 
Z(K)=Y (K)*H( I-N02, J,K) 

23 ZCK+NH02) =Y (K+NH02)»H( 1-N02, J,K} 

2( 1 1=Y{ I )*H{ I-N02,J, 1 ) 

ZINH21 1 =Y (NH21 ) »H { I -N02, J,NH21 1 
GO TO 21 

2a DO 25 K=2,NH02 

2(K)=Y(K:*H(I, J-N02 , K > 

25 Z(K+NH02 J=Y<K+NH02 )*H( I , J-N02,K ) 

Z< I ) =Y( I )»H( I , J-N02. ’ ) 

2(NH21 )=Y(Nh 21 )»Hfl ,J-N02,NH21 ) 

GO TO 21 

26 00 27 K=2. NH02 
2{K)=Y(K]*H{ I-N02.J-n02,K) 

27 2(K+NH02 l=YlK+NH02l*H(I-N02,J-N02,)C) 
Z t 1 ) =Y ( 1 1 »H ( I -N02 , J-N02 , 1 1 

2(NH21 ) =Y (NH21 ) *H ( I -N02, J-N02,NH2l 1 
21 CONTINUE 

call GETSET C a, !3A ) 

CALL FTPANS (A, I3A ) 

DO 23 K=1,nH02 

28 RHOf I , J.K l=Yti<) 

20 CONTINUE 

CALL GETSET ta, I2A ) 

DO 29 K=i,NH02 
DO 29 J=1 ,N 
DO 30 I = 1 ,N 

30 Zf I )=RHO{ I , J.K) 
call ^'TRANSia, I2A) 

DO 31 1=1 .N 

31 RFO( I , J.K )=Yt I ) 

29 CONTINUE 

DO 32 K=1 ,NH02 
DO 32 1=1 ,n02 
' DO 33 J=1 .N 

33 Z t J)=RHO( I . J.K ) 

CALL FTRANS(a,12A) 

DO 34 J=I.N02 

34 RHOf I . J.K)=Y t J) 

32 CONTINUE 
RETURN 
END 
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■ TABLE A2 

OVERLAYS FOR CALCULATING THE THREE-DIMENSIONAL GRAVITATIONAL 
POTENTIAL USING CORE AND DISK STORAGE 
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n o n n 


page 6 

^ QUAIJTV 


c the following is thl section of an- initializing overlay in which constants 

C (RELATED to Thfc DIMENSIONS OF THE PHI ( OENS I TY/POTENT 1 AL ) ARRAY ARE CO''- 
C PUTEO. IT IS CALLED ONCE AT THE BEGINNING OF A PROGRAM run. IN THIS ^ 

C LISTING THE VALUES OF 12A. ISA AND THE DIMENSION AND LABELED COMMON 
C STATEMENTS ARE SET FOR AN ACTIVE PHI ARRAY DIMENSIONED 6A BY 6A BY 16. 

I2A = 7 

I3A = 5 

I2B= 12A-1 

I3B=I3A-1 

N=2**I2A 

N02=N/'2 

N21=N02-H 

NOAsN/A 

N3A=NO2•^N0A 

NH=2**I3A 

NH02=NH/2 

NH2I =NH02-i-l 


c ****-*-»-»*****-*-*R'*»R-*ir******»*«-********;<.**A.*.i(.**j.j(.* + ***.:*j.*^^^*j.*^^.*;fj^jj.^jtj.jt*^^ 

c the following IS The section of the star advancing overlay in which cnuN<s 

C of the PHI ARRAY (CONTAINING THE DENSITY MESH) ARE -.RITTEN ONTO D1S< FILES 
C 1.2i5 AND 6. THt star advancing OVERLAY IS CALLED ONCE PER TIME STEP. 
DIMENSION PHI (6A.6A.16).OI 132.32.16) 

DO 520 K=!.nH 02 
DO S20 U=1.N0A 
DO 520 1=1, NO A 
520 01 < I , J,K )=PHI t I . J.K ) 

WRITEU > 01 

rewind I 

DO 525 K=1.NH02 
DO 52S J=1 .NOA 
DO 525 1=1 ,NOA 

525 01 I I . J.K)=PhI I I .NOA-t-J.K) 

WRiTEtS) 01 

rewind 5 

DO 330 K=I .NH02 
DO 530 J=I ,NOA 
DO 530 1=1, NOA 
530 01 ( I , J.K )=PHI (N04+! . J.K) 

WRITE<2) 01 
REWIND 2 
DO 535 K=1,NH02 
DO 535 J=1.N04 
DO 535 1 =1 ,N0A 

535 01 <I . J,K)=PHI (NOA-I-I ,NOA+J,<) 

WRITE<6) OI 
REWIND 6 
C 
C 
C 

c the following IS The section of the star advancing overlay in which chun<s 
C of the PHI array (CONTAINING THE POTENTIAL MESH) ARE READ FROM DISK FILES 
C 1,2,5 And 6. 


001 

002 

003 

OOA 

005 

006 

007 

008 

009 

010 
Ol 1 
012 
013 
Ol A 

015 

016 
017 
OIS 

019 

020 
021 
022 
023 
02a 

025 

026 

027 

028 

029 

030 

031 

032 

033 
03a 

035 

036 

037 

038 

039 
OAO 
OAl 
0A2 
0A3 
OAA 
0A5 

046 

047 

048 

049 

050 

051 

052 

053 
05A 
055 
066 
037 
058 
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. dimension phi (64,64. IS) ,01 (32.32, 16 ) 
read ( 1 ) 01 
rewind 1 

DO 30 K=1 ,nh02 
DO 30 J=1 ,n04 
DO 30 1=1 ,n04 
30 PHK I . J,K )=0I { I ,J,K) 

R£A0(5) 01 
REWIND 5 
DO 40 K=i ,NH02 
DO 40 3=1 ,n04 ■ 

DO 40 1=1 ,no4 

40 PHI ( I ,N04+3,K )=0I ( I . J,K ) 
read (2) 01 
REW I NO 2 


DO 50 K=1,NH02 
00 50 3=1 ,n04 
DO 50 1 = 1 ,n04 

30 PHI (N04+I ,3,K )=Ol ( I ,3,K) 
READ{6) OI 
rewind 6 

DO 60 K=l ,NH02 ' 

DO 60 3=1 ,N04 
DO 60 1=1 ,n04 

60 PHHN04+I ,N04+3,K 1=01 { I ,3,J<) 
C 
c 
c 


following is The geth overlay, which computes i'.o stores the trans- 

C -ORMtD GREENS FUNCTION. IT IS CALLED ONCE AT Ti== cEGINNIng'of A PROGRAM 
C RUN • 


OVERLAY! IFILE.4.0) 

PROGRAM GETH 

c This overlay performs a cosine analysis of the three-dimensional greens 
C function array, it then writes chunks of this array on disk file 9 IN the 

C OHDER IN WHICH THEY WILL BE READ INTO THE HH ARRA'' DsjPIng THE GETPH I 
C OVERLAY. VALUES FOR l=N/a+! AND J=N/2+l ARE TRANSFERRED TO THE hn21 ARRAY 
C WHICH IS IN COMMON WITH THE GETPH! OVERLAY, 

COMMON/ ALLCOM/n »N02 «N2 1 .N04 .N34.NH.NH02 .NH2 L • 1 2A . I2S » I3A . I3B 
C0MM0N/HN21C0M/HN2I (65.17) 

COMMON 2(1025). Y ( I 025 ) 

O I MENS ION H ( 65 , 65 . 1 7 ) . HH ( 32 . 32 , 1 7 ) 

RNI = 1 ,/(N*-N*NH) 

DO 1 K=I,NH21 
DO 1 3= 1 , N2 1 
DO 1 1=1, N2 I 

Rl = (K-l )* (K-1 ) + (3-l )* (3-1 l + ( I-l >♦ ( I-l ) 

IF(RI.LT*1. 1 Rl = j p 
H( I .3.K )=RN!/SQ«T(RI ) 

1 CONTINUE 

call GETSET(2. I2S ) 

DO 2 K=1 ,NH21 
00 2 3=1 ,n21 
DO 3 1=1 ,N21 
3 2( I )=H( I ,3,K) 


059 

060 
061 
062 

063 

064 

065 

066 

067 

068 
069 
OTO 

071 

072 

073 

074 

075 

076 

077 

078 

079 

080 
081 
032 

083 

084 

085 

086 

087 

088 

089 

090 

091 

092 

093 

094 
093 

096 

097 
09B 

099 

100 
101 
1 02 
103 
t-04 

105 

106 

107 

108 
1 09 
1 1 0 
1 I I 
112 
1 13 
1 14 
I 15 





CAUL FTRANS (2 .'128 ) 

DO 4 I=1.N21 

4 H( I , J,K)=Y( I t 

2 CONTINUE 

DO 5 K=2 .NH21 
DO 5 1=1 .N21 
DO 6 J=1 .N21 

6 2( J)=H( I , J.K) 

CALL FTHANSI2, [2B ) 

DO 7 J=1 »N21 

7 H( r . J.K J=Y( J 1 

5 CONTINUE 

CALL GETSET(2, I2B ) 

DO 10 J=I ,n2! 

DO 10 I = I , N2 ! 

DO e K=1 ,NH2! 

a z (K)=H ( I ,j.K) 

CALL FTRANS(2* I3BJ 
DO 9 K=1 .NH21 

9 HU.J*K)=Y(KJ 
10 CONTINUE 

DO 30 1=1 ,N04 
DO 30 J=1 »N04 
DO 30 K=1 ,NH21 
30 HHI I , J.K )=H( I • JtK ) 
WR1TE19) HH 
DO 35 1=1 .N04 
DO 35 J=1 ,N04 
DO 35 K=1 ,NH21 
35 HH( I . J.K1=H( I .N04+J.K) 
WRITEI9J HH 
DO 40 1=1 .N04 
DO 40 J=I ,N04 

DO 40 K = 1,NH21' 

40 HH(I,J.K>=H< I.J.K) 

WRITE(9) HK 
00 45 I=l ,N04 
DO 45 J=1 .N04 
DO 45 K=1 ,Nri21 
45 HH(I.J.K)=H{N04+t .J.K) 
WRITE<9) HH 
DO SO 1=1 .n04 
DO 50 J=1 ,N04 
DO 50 K=1 .NH21 

50 HH t I . J.K)=H (N04+I .N04+J.K 1 
WR1TE<9 ) ,HH 
DO 55 1=1 .N04 
DO 55 J=1 .N04 
DO 55 K=1 .NH2I 
55 HH ( I . J.K ) =H tN04+I . J »K ) 
WRITEt9) HH 
rewind 9 
DO 15 K=1 ,NH21 
DO IS I=l .N21 
IS HN21 ( I .K)=H( I .N21 .K) 

RETURN 

END 
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c 

c 

c 

c*-*^»-»******»*** *•■»•»***»*)>♦■*■***»*»■»*»■♦*»•*■»•***♦*■**********. ■»■*♦■*.***■*■*♦»■**»■»•**■»»» 
c the following is the getphi overlay, which computes the potential fesh. 

C IT REPLACES CHUNKS OF DENSITY STORED ON DISK FILES 1«2»S AND 6 WITH 
C CORRESPONDING CHUNKS OF THE POTENTIAL MESH. IT IS CALLED ONCE PER TIME 
C STEP. 


overlay ( GF I LE . S , 0 } 

PROGRAM GETPHI 

C Ti-<is overlay solves for the potential mesh {dimensioned n/2 by N/2 by 

C NH/2 > DUE to A DENSITY MESH (DIMENSIONED N/2 BY NKE BY NH/2 ) BY DOING A 

C PERIODIC ANALYSIS OF THE DENSITY AND THEN A PERIODIC SYNTHESIS OF THE 
C PRODUCT OF THE TRANSFORMED GREENS FUNCTION (DIMENSIONED (N/2+1 ) SY (N/2+I ) 
C BY (NH/'2+I)) AND The transformed DENSITY. FORMALLY SPEA<ING, EACH OF THE 
C transforms (EXCEPT THE COSINE ANALYSIS OF THE GREENS FUNCTION, WHICH IS 
C performed in THE GETH OVERLAY) REQUIRES AN ARRAY DIMENSIONED N BY N BY 
C NH. TO REDUCE CORE STORAGE THIS OVERLAY PERFORMS THESE TRANSFORMS IN 
C CHUNKS BY THE alignment OF FOUR SMALLER ARRAYS NAMED RHOI , HH02 , RH03, AND 
C PHOA, EACH OF WHICH IS DIMENSIONED N/4 BY N/4 BY NH/2. THE CHUNKS OF THE 
C lower half (I .LE. 2 .LE. NH/2) OF THE EXTENDED ARRAY NOT IN CORE AT ANY 
C ONE TIME are STORED ON DISK FILES 1 THROUGH 0. THE FOLLOWING ARE TWO TOP 

C VI-wS Or THE LOWER HALF OF THE EXTENDED ARRAY. BOTH OF THESE VIEWS 

C DESIGNATE the CHUNKS AS IROW AND JCOLUKN. IROW I AND 2 OF JCDLUMN 
C 1 AND 2 CONSTITUTE THE ACTIVE MESH. In THE DIAGRAM ON THE LEFT THE 
C NUMBERS WITHIN THE CHUNKS OF JCOLUMN I AND 2 INDICATE THE DISK FILES ON 
C which THOSE CHUNKS ARE STORED. (NO DISK STORAGE IS REQUIRED FOR JCOLUMN 3 
C OR 4. ) REFERRING TO THE DIAGRAM ON THE RIGHT, THE NUMBERS WITHIN THE 
C CHUNKS are the ORDER IN WHICH CHUNKS OF THE TRANSFORMED DENSITY ARE 


MUTIPLIEO (ELEMENT BY ELEMENT) BY THE APPROPRIATE PORTION OF THE 

transformed greens function which has been read from disk file 9 Into 
array HH(N/4,N/4,NH/2+1 ). (AN EXCEPTION IS THE SET OF TRANSFORMED GREENS 
FUNCTION BOUNDARY VALUES FOR I=N/2+I AND J=N/2+I WHICH REMAIN AT ALL TIMES 
IN COMMON IN THE array HnE 1 (N/2+ 1 ,NH/2+I ) . ) A PLUS IN A CHUNK INDICATES 
THAT NEW VALUES MUST BE READ INTO ARRAY HH BEFORE THAT CHUNK IS MULTIPLIED 
BY hH. THIS SYSTEM MINIMIZES PERIPHERAL PROCESS TIME SY UTILIZING THE 
periodicity of THE TRANSFORMED GREENS FUNCTION. 


TWO 


TOP VIEWS OF LOWER HALF OF EXTENDED MESH(N BY N SY NH/2) - IROW I 
AND 2 OF UCOlUMN 1 AND Z CONSTITUTE THE ACTIVE MESH(N/2 BY N/2 BY 
NH/2). THE directions are XU) AND OMEGAX ( I ) - DOWN ON PAGE, 

Y(J) and OMEGAY(J) - TO RIGHT ON PAGE, Z(K) AND OhEGAZ(K) - OUT OF 
PAGE. 


JCOLUMN 
12 3 4 

**♦*»*»*■)(•■»•*****■»* 

IR0W=1 * I » 5 * * 

***♦■*• 
IR0W=2 * 2 *■ 6 * * * 

^ * 
IR0W=3 *3*7* *■ * 

♦ » * ♦ * 

1R0W=4 * 4 * 8 ♦ * 

disk: files on which chunks 
are stored 


JCOLUMN 
1 2 3 4. 

« 4. if. .f, « ^ ^ 

IR0W=1 * 3 * Z ^ ^ * 

♦ * * 

IR0W = 2 ♦ 9 *-11 -*10 *12 ■*" 

* .f. » «■ # 

IR0W=3 ♦7*5*8*6* 

-»« if’ 4- if- ^ * 

^ * * * * 

IR0W=4 +15 *13 *16 *14 * 

ORDER IN WHICH CHUNKS ARE 
MULTIPLIED BY APPROPRIATE 
PORTION O- transformed 
GREENS FUNCTION 


C 

c 


173 

174 

175 
1 76 

177 

178 

179 

lao 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 
2C3 

204 

205 

206 

207 

208 

209 

210 
21 1 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 
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WKiiiiNALf lo 

OE POOR QUALITY 


COMMCN/ALLCOM/N.NOa.Nai .N04 .N34,NH.NH02,NH2I . 12 A, I SB • I3A. 13B 
COMMON/TRANCOM/RhOI <32.32. 16) .RH02 (32.32. 16 ) . RH03 ( 32 .32, 16 ) . 

I RH04<32.32,I6).HH(32.32.17) oc.ioj. 

CO(.'MON/'HN21COM/.riN21 (65.17) 

thc initializing overlay or star advancing overlay stores the density 

‘ ^ ^COLOMN=I ON disk files I AND 2 RESPECTIVELY 

RPO, 6 respectively, the GETPHI OV^Ri_AY 

Pn-^WT-^ density ON these DISK FILES WITH THE CORRESPONDING VALU=-S OF 
potential which are then used in the star advancing overlay. thIs li 
accomplished through calling subroutines anlxi jcolumn) . anlsyncirow) and 
SYNXI JCOLUMN) AS detailed BELOW. 


subroutine ANLX(JCOLUMN) reads respectively IROW I AND 2 from the 

Plues - t and 2 FOR UCOLUMN=t, - 5 AND 6 FOR JC0LUMN=2. 

I PERrORMS A PERIODIC ANALYSIS IN THE X DIRECTION OVER JCOLUMN FOR 

1-1, N AND WRITES THE RESULTS RESPECTIVELY FOR IROW 1.2.3, AND 4 ON THE 

UC0lCmn= 2.'’‘^" UCOLUMN = l, - 5.6.7, AND 8 FOR 

CALL ANLX ( 1 ) 

CALL ANLX (2) 

C SI^ROUTINE ANLSYNCIROW) READS RESPECTIVELY JCOLUMN 1 AND 2 FROM THE 

7 FOR^'Rnw‘'3^^ - t and S FOR 1R0W=I. - 2 AND 6 FOR IRoS^sI - 3 AND 

** ^ ® IROW=4. IT TH£N PSRFOR.'-S A PEftlOOlC ANALYSIS 

^ DIRECTION OVER IROW FOR J=1.N. FOR EACH CHUNK I ^THEN PERFORmI ^ 
PcRtODIC ANALYSIS IN THE Z DIRECTION FOR K=I,NH, ELEMENT BY ELEMENT 

with a similarly SHAPED CHUNK OF THE TRANSFORMED GREENS 
RFSn, A PERIODIC SYNTHESIS IN THE 2 DIRECTION FOR K = i,NH. THE 

RESULT for K=!,NHXS IS THEN PERIODICALLY SYNTHESIZED IN THE Y DIRECTION 
J = 1.N. THIS last result for JCOLUMN I AND 2 IS WRITTEN 

FOLLOWING DISK FILES - I AnO 3 FOR IROW=I. - 2 AND 6 
FOR IR0W=2. - 3 AND 7 FOR IROW=3. - 4 AND 8 FOR IROW= 4 . THE ORDER IN 

called for IROW 1 THROUGH 4 MINIMIZES READING FROM DISK 
OF THc transformed GREENS FUNCTION AS MENTIONED ABOVE, 

call anlsyn cl) 

CALL anlsyn (3) 

CALL ANLSYN (2) 
call ANLSYN (4) 

SUBROUTINE SYNX( JCOLUMN) READS RESPECTIVELY IROW 1.2.3. And 4 FROM THE 
PILES - 1,2,3, ANO 4 FOR JCOLUMN=l. - 5,6,7, ANO 8 FOR 

PERFORMS A PERIODIC SYNTHESIS IN THE X DIRECTION OVER 
COLUMN FOR J=1,N, IT THEN WRITES THE RESULT RESPECTIVELY FOR IROW 1 ANO 

- 1 and 2 FOR JC0LUMN=i, - 3 °nO 6 FOR 

CALL SYNX( I ) 
call SYNXC2) 

RETURN 

end 

SUBROUTINE ANLX ( JCOLUMN ) 

C0MM0NZALLC0M/N.N02.N2I .N04.N34.NH.NH02.NH2l , I2A, 120, J3A, I3B 
COMMONZTRANCOM/RhOI ( 32 . 32 . l 6 ) .RH 02 C 32 . 32 , 1 6 ) . RH03 (32.32,16). 

1 RH04 ( 32, 32, I 6 ) , HH ( 32. 32 , 17) 

COMMON 2(1023), Y(1025) 

IF (JCOLUMN. EO. 2 ) GO TO 2 
RE AO (I ) RHOI 
REW I ND 1 
read (2) RH02 

rewind 2 
GO TO 3 

2 CONTINUE 
read (5) RHOI 

rewind 5 
Read ( 6 ) rho 2 
rewind 6 

3 CONTINUE 

CALL GETSET(3. I2A) 

DO 10 K=! ,NH02 
DO 10 J=1 ,n04 
OO 51=1 ,N04 
Z( I )=RH01 ( I . J.K) 

Z(N04+i )=RH02(1 .J.K) 


244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 

266 
267 
260 
269 
370 

271 

272 

273 

274 

275 

276 

277 

278 

279 

280 
281 
282 

283 

284 

285 

286 
237 
288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 
3 1 4 
315 
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2CN02+I )=0. 

5 Z(N34+I )=0. 

call FTRANS (3i I2A ) 

DO 1 a 1 = 1 , N‘04 
RHOI (.! , J,K) = Y«J ) 

RH02C I . J,K)=y tNOA+I 1 
RH03t I . J.K)=Y(N02+I > 

10 RHOA ( I , J.K) =Y (N34+I ) 

IP'fJCOLUMN.EQ.a) GO TO 12 
WRiTil I ) RHOI 

R£wmo 1 
WRITE (2) RH02 
Rewind 2 
WRITS<3) RH03 
rewind 3 
WRITE! A) RHOA 
REWIND A 
GO TO 15 
12 CONTINUE 

WRITE (5) RHOI 

Rewind s 
WRIT e<6) RH02 
rewind 6 
WRITEI7) RH03 
rewind 7 
WRITEtSl RHOA 

rewind 8 

15 RETURN 
END 

subroutine ANLSYN(IROW) 

COMN0N/ALLC0M/N.NO2.N21 .NOA.NOA.NH, NH02 ,NH2l . I2A. 123, I3A, I3S 
COMH0NyTRANC0MxRH0l<32,32.i6>.RHO2(32.32.ir,.RHO3 32 32^ 6)! 
I RH0A<32.32,16),HH(32.32.17) ‘ =>2 ,32. 1 6 ) . 

COMMONZHN2ICOM/HN21 (65,17) 

COMMON Z(1025), Y(102S) 

GO Toil, 2, 3, A) I ROW 

1 CONTINUE 

REAO(I) RHOI 

rewind I 
read (3) RH02 
rewind 5 
GO TO S 

2 CONTINUE 
Read (2) rhoi 

REWIND 2 
read (6) RH02 
rewind 6 

GO TO 5 

3 Continue 
Read (3) rhoi 
rewind 3 
read (7 J RH02 
rewind 7 

GO to 5 
A CONTINUE 

R£ad(A) Rhoi 
Rewind a 
Read ( 3 j rho2 
rewind 8 
5 continue 

CALL GETSET(3. I2A > 

DO 10 K=l ,NH02 
DO 10 1=1, NOA 
DO 7 J=1 ,N04 
ZC J)=Rh 01 ( I , J,K ) 
z ( NOA+ J ) = rho 2 r I , j ,k ) 

Z(N02+J)=0. 

7 Z(N3A+J)=0. 

CALL FTRANSO. I2a ) 
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UKltrUNAL I'AGJU Jis 

OF POOR QUALrry 


DO 10 J = ! .N04 _ -jas 

RHOl ( I . J.K)=Y( J)‘ 386 

RH02( I . J»K)=YCN04+J) 387 

RH03f I . J«K)=Y(N02+J) 388 

O RH04 ( I . JtK )=Y (N34+J > 389 

GO TO<30.49,75.75 ) I ROW 390 

CONTINUE 391 

CONTINUE 392 

■READ (9) HH 393 

iO JC0LUMN=1 394 

DO 70 1=1 .N04 . 395 

DO 70 J=1 ,N04 396 

DO 52 K=1.NH02 397 

ZtK)=RH0! ( 1 , J.K ) 398 

i2 2INH02+K)=0, 399 

CALL GETSET(3, I3A ) 400 

call FTRANS <3. 13A 1 401 

IF t IR0W.NE.3 ) GO TO 300 402 

IFCl.NE.llGO TO 300 403 

LL=J 404 

GO TO 200 405 

>4 DO 70 K=I .NH02 406 

’O RHOl t i • J»K}=YCK ) 407 

GO TO 100 408 

■4 continue. 409 

read (9) HH 410 

’5 JC0lUMN=2 411 

DO 95 I =l .N04 412 

DO 95 J=1 .N04 413 

DO 77 K=l.NH02 414 

2(K).=RH02 f I ) 415 

'7 2(NHO2+K)=0. 416 

CALL GETSETO. 13A ) 4 17 

CALL FTRANSOt I3A J 418 - 

IF(IROW.NE.3j GO TO 300 419 

IFII.NE.llGO TO 300 420 

LL=N04+J 421 

GO TO 200 422 

■g DO 95 K = 1,NH02 423 

5 RH02I I « J.K>=Y(K> 424 

GO TO 12S 425 

0 JC0LUMN=3 426 

DO 120 1=1, N04 427 

DO 120 J=1,N04 428 

DO 101 K=1,NH02 429 

2 (K)=HH03 t 1 430 

'1 2<NH02+K)=0. 431 

CALL GETSETO, I 3A ) 432 

CALL FTRANSO, OA ) 433 

GO 70(103,105,107,115) I ROW 434 

3 IF(J.NE,l)GO TO 300 435 

LL=I 436 

GO TO 200 437 

5 IF(J,NE.l )GO TO 300 438 

LL=N04+I 439 

GO TO 200 440 

■7 IF ( 1 ,NE.l .A nd. J.nE. l)GO TO 300 441 

IF ( I ,EQ,l .AND.J.EQ. 1 )CO TO 111 442 

lF(I,EQ.I)GO TO 109 443 

LL=I 444 

GO TO 200 445 

9 LL=J 446 

Go TO 200 447 

1 LL=N21 448 

GO TO 200 449 

5 IF(J.NE.l) GO TO 300 450 

LL=N04+I 451 

GO TO 200 452 

7 DO 120 K=1,NH02 4S3 

:0 RH03( I , J,i<)=Y{K ) 454 

GO 70(74.74.400,390) I ROW 455 

;5 JC0LUMN=4 456 



ORIGINAL Page is 
OF POOR QUALITY 


DO 145 1=1, N04 
DO 145 J=1,N04 
DO 137 K=1.NH02 
Z(K)=RH04(I ,J.K) 

137 Z(NH03+K)=0. 

call GETSETO, I3A ) 

CALL FTRANS{3.I3A) 

IF ( IROW .NE.3 ) GO TO 300 
IF( I .NE.l )G0 TO 300 
LL=N04+J 
GO TO 300 

139 DO 145 K=1.NH02 
145 RH04( I , J,K)=Y(<) 

GO TO (400,400,49,49) IROW 
300 DO 305 K=3.NH03 

Z(K)=Y(K)*HN21 (LL,<J 
305 2<NH03+K ) =Y ( NH03+K ) 4HN31 (LL,K) 

Z(l )=Y(1 )»HN21 (LL,1 I 
Z (NH21 >=Y (NH31 )*HN31 (LL,NH31 ) 

GO TO 310 

300 DO 305 K=3,NH03 

2<K)=Y<K)*hh< I ,3,K) 

305 2 tNH02+K ) =Y ( NH03+K ) *HH C I , J,K ) 

Ztl )=Y ( I )*HH( 1,3, 1 ) 

Z(NH21 )=Y(NH31 )*HH( I .J.NH31 ) 

310 call GETS£T{4,I3A] 
call FTRAnS(4. I3A 1 
GO TO ( 54,79, 1 17, 139) JCOLUMN 
390 REWIND 9 
400 CALL GETSET(4,I3ai 
DO 410 X:=I.NH03 
DO 410 1=1 ,N04 
DO 405 J=I,N04 
2( J)=RH01 ( 1 , J.K ) 

Z(N04+J)=RH02( I ,J,K) 

2 (N03+J J =RH03 ( 1 ,K ) 

405 ZtN34+J)=RH04tI ,J.K) 

CALL FTRANS{4, I3A ) 

DO 410 J=1 ,N04 
RHOl { I , J.K»=Y( J) 

410 RH03 [ I , J,K) = Y1N04+J ) 

GO TO<415,430,43S.430 ) IROW 
415 CONTINUE 

write n ) RHOl 
rewind I 
WH!TE(5) RH03 
REWIND 5 
GO TO 435 
420 CONTINUE 

WRITE (2) RHOl 
REWIND 2 
WRITE (6) RH02 
rewind 6 

GO TO 435 
425 CONTINUE 

WRITE(3) RHOl 
REWIND 3 
WRITEI7) RH02 
REWIND 7 
GO TO 435 
430 CONTINUE 

WR1TEI4) RHOl 

rewind 4 

WRITE (8) RH02 
REWIND 8 
435 RETURN 
ENO 

subroutine SYNX( jcOLUMN) 

COMMONS ALLCOM/N ,N02 ,n21 ,N04 , N34 ,NH ,NH02 ,NH2 1 ,I2A,I2B, I3A, I 38 
COMHON/TR ANCOM/RhOI (32,32,16) , RH02 (32,32,16) ,RH03 (32,32,16), 
1 RH04 (32,32, 16 ) ,HH(32,32 , 17 ) 

COMMON Z( 1025) ,Y( 1025) 

IF(UCOLUFN.EO,2) GO TO 1 


457 

458 
439 
4*60 

461 

462 

463 

464 

465 

466 

467 

468 

469 

470 

471 

472 

473 

474 

475 

476 

477 

478 

479 

480 

481 

482 

483 

484 
^485 

486 

487 
483 

489 

490 
49! 

492 

493 

494 

495 

496 

497 

498 

499 

500 

501 

502 

503 

504 

505 

506 

507 

508 

509 

510 
5! 1 

512 

513 

514 

515 

516 

517 

518 

519 

520 
52! 

522 

523 

524 

525 

526 

527 

528 

529 
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ORIGINAL PAGE Ib 


REAOJl) RHOt 

OE POOR QUALITY 

530 

REWIND 1 


531 

read (2) RH02 


532 

REWIND 2 


533 

READ(3) RH03 


534 

REWIND 3 


535 

REAOCA) RH04 


536 

REWIND- A 


537 

GO TO 2 


538 

1 CONTINUE 


539 

read ( 5) RHOl 


540- 

REWIND 5 


541 

read 16 I RH02 


542 

REWIND 6 


543 

READ(T> RH03 


544 

REWIND 7 


545 

read (8) RH04 


546 

REWIND e 


547 

2 CONTINUE 


543 

4 CALL GET5ETI4, 12A ) 


549 

DO 10 K=1 ,NH02 


550 

DO 1 0 J=1 ,N04 


55! 

DO 5 I=1.N04 


552 

2(1 )=RH01 ( I .J.K) 


SS3 

Z(N04-H >=RH02< I .J .K) 


554 

ZCN02-H )=RH03( I tJtK) 


555 

5 ZCN34-(-I !=RH04{ I .J.K) 


556 

CALL FTRANS (4. I2A ) 


557 

DO 10 1=1 .N04 


55S 

RHOI ( I . J.K)=Y( I ) 


559 

10 RH02 ( I . J.K)=Y(N04-H ) 


560 

IF (JCOLUMN.EO.2 ) GO TO 12 


561 

WRITEU) RHOl 


562 

REWIND ! 


563 

WRITE(2) RH02 


564 

REWIND 2 


565 

GO TO 15 


566 

12 CONTINUE 


567 

WRITEO) RHOl 


568 

REWIND S 


569 

WRITEI6) RH02 


570 

REWIND 6 


571 

15 RETURN 


572 

END 
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TABLE A3 


Array Dimensions 
(Program of Table A2) 


Array name 

General 

dimensions 

Dimensions used . 
in actual runs 
and listing of 
Table A2 

1 

Dimensions used 
in Figs. A1-A6 

Overlays in which 

dimensioned 


(note 1) 

Star adv. 6ETH 
and initl . 

GETPHI 

PHI (active) 

n X n X h 

64 X 64 X 16 

8 X 

8x4 

x 


01 

(n/2) X (n/2) x h 

32 x 32 X 16 

4 X 

4x4 

X 


H 

(n+1 ) X (n+1 ) x (h+1 ) 

65 X 65 X 17 

9 X 

9x5 

X 


HH 

(n/2) X (n/2) x (h+1) 

32 X 32 X 17 

4 X 

4x5 

X 

X 

HN21 

(note 2) 

(n+1 ) X (h+1 ) 

65 X 17 

9 

x 5 

X 

X 

RH01,RH02 

RH03,RH04 

(n/2) x (n/2) x h 

32 X 32 X 16 

4 X 

4x4 


X 

Extended 

PHI 

(2n) X (2n) x (2h) 

128 X 128 X 32 

16 X 

16 X 8 

not actually dimensioned 
(note 3) 


Cn 


Note 1: The notation a x b x c represents the array dimensions of the subscripts x, y and z, 

respectively, (or the subscripts c, n and respectively, of the transformed array) such that 
a X b X c equals the total number of array elements. The Fortran variables N and NH are equal 
to 2n and 2h, respectively. 


*\j 

H. 


elements. Its first 


Note 2: HN21 is a two-dimensional array containing a boundary plane of .. ^ ^ 

subscript corresponds to 5 or n equivalently, while its second subscript ’corresponds to 5 . 


Note 3: While the program uses smaller arrays in order to avoid dimensioning the (2n) x (2n) x (2h) 
extended PHI array of Fig, 1, its mathematical existence is necessary for the Fourier solution of the 
potential of an isolated galaxy. 


r^TXT A T p&aR TS 


ORIGINAL PAGE Ib 
OB PCX)R QUALITY 

TABLE AA 


*\j 

Storage of the Fourier Transformed Green's Function H 

on Desk File 9 


{Program of Table- A2) 


Record No. Storage sequence within Use sequence within GETPHI 

of file 9 GETH- overlay (Note 1) ’ overlay (Note 2) 


1 

2 

3 

4 

5 

6 


A 

B 

A 

C 

D 

C 


(1J),(1,3) 

(1,2), (1.4), (3, 2), (3,4) 

(3.1) , (3, 3) 

( 2 . 1 ) , ( 2 , 3 ) 

(2.2) , (2,4), (4, 2), (4,4) 
(4,1), (4, 3) 


Note 1: Within the GETH overlay, this is the location in the H array (as 

designated by letters A-D of Fig. A4) from which "do loop" transfer is made 
to the HH array followed by writing on the Indicated record of disk file 9. 

Note 2: Following reading of the indicated record of disk file 9 into the 

HH array within the GETPHI overlay, this is the sequence of locations in 
the extended PHI array (as designated by "chunks" (IROW.JCOLUHN) of Fig. A2) 
upon which z-direction one-dimensional^array operations are performed. These 
operations include multiplication by H, the appropriate portion of which is 
now contained in the HH array. This method minimizes reading of file 9 by 
using the periodicity and symmetry of ?f. 


52 



TABLE AS 


Outline of the GETPHI Overlay 
(Program of Table A2) 


(Refer to Fig. A6 for orientation of arrays RHOl, RH02, RH03 and RH04 
and to Fig." A2 for file, numbers corresponding to the "locations" of 
these arrays.) 

Listing line Nos. 
of Table A2 


A. CALLANLX(l): Fig. A6(a). 

1. Read files 1 and 2 into RHOl and RH02, respecti vely. 299-302 

2. Set RH03=RH04=0. 316-317 

3. Perform Fourier transform in x-direction over RHOl, 310-323 
RH02, RH03 and RH04: p -> 

4. Write RHOl, RH02, RH03 and*RH04 onto files 1, 2, 3 325-332 

and 4, respectively. 


B. CALL ANLX(2): Fig. A6(b). 

1. Read files 5 and 6 into RHOl and RH02, respectively 

2. Same as steps A. 2 and A. 3 

3. Write RHOl, RH02, RH03 and RH04 onto files 5, 6, 7 
and 8, respectively. 


305-308 

335-342 


C. CALL AMLSYN(l): Fig, A6(c). 


1 . 

2 . 

3. 

4. 

5. 


Read files 1 and 5 into RHOl and RH02, respectively 
Set RH03=RH04=0 

Perform Fourier transform in y~direction over RHOl, 

RH02, RH03 and RH04: p^ 

^>y » z n ,z 

Read record 1 of file 9 into HH 
For each one- dimensional array in z-direction of 
which RHOl is composed: 
a. Transfer to one-dimensional array Z, 
dimensioned at least 2h+l 
Set Z=0 for z >. h ' ’ 

Perform Fourier transform in z-direction over 


b. 

c. 


Z for 0 ±z ± 2h-l 
-in one-dimensional 

Y by ft. 


with the result appearing 
array Y: 


d. Multiply 


353-356 

382-383 

376-389 

393 


ft. 


e. 


& 

^ to form = ft 

Perform inverse Fourier transform in z-direction 
over Y and store result for 0 <. z <, h - 1 
in RHOl: ft. „ ^ ft 






6. 

Repete step C.5 for RH03 

426-454,471-483 

7. 

Read record 2 of- file 9 into HH 

410 

8. 

Repete step C.5 for RH02 and RH04 

41 1 -424 , 456-469 , 477-483 

9. 

Perform inverse Fourier transform in y-direction 
over RHOl, RH02, RH03 and RH04: „ 

486-497 

10. 

Write RHOl and RH02 onto files 1 and 5, 
respectively. 

500-503 
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Listing Line Nos. 
of Table A2 

D. CALL ANLSYf|(3): Fig. A'6(e). 

1. Read files 3 and 7 into RHOl and RH02, respectively. 

2. Same as steps C.2-C.9 except for sequencing of 
reading tape 9 into HH and the z-di recti onal 
operations. Table A4 details this sequencing. 

3. Write RHOl and RH02 onto files 3 and 7, respectively. 

E. CALL ANLSYN(2): Fig. A6(d). 

1. Same as step D except that files 2 and 6 correspond 
to RHOl and RH02, respectively, for read and write 
operations. 

F. CALL ANLSYN{4): Fig. A6(f). 

1. Same as step D except that files 4 and' 8 correspond 
to RHOl and RH02, respectively, for read and write 
operations. 

G. CALL SYNX(l): Fig. A6(a). 

1. Read ’files 1,2, 3 and 4 into RHOl, RH02, RH03 and 530-537 

RH04, respectively. 

2. Perform inverse Fourier transform in x-direction 550-560 

over RHOl, RH02, RH03 and RH04: y ^ y z 

3. Write RHOl and RH02 onto files 1 and 2, respectively. 562-565 

H. CALL SYNX(2): Fig. A6(b). 

1. Read files 5, 6, 7 and 8 into RHOl, RH02, RH03 and 540-547 

RH04, respectively. 

2. Same as step G.2 

3. Write RHOl and RH02 onto files 5 and 6, respectively, 568-571 


365-368 

512-515 
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Z,(?) 



Figure Al.- PHI array (active), which contains the galactic density/potential mesh, and 
the extended PHI array, which is required for the Fourier potential solution of an 
isolated galaxy. Each x-, y-, or z-axis represents the following: (a) the x-, y-, or 
z-spacial direction; (b) the unt runs formed array subscriipt x, y, or z; and (c) the x-, 
y-, or z-diroction transformed array subscript q or i,, respectively. For clarity 
in this and the following figures, the PHI array is dimensioned n x n x h «= 8 x 8 x d 
^ while in the program as listed in Table A2 and as actually run it is dimensioned 

64 X 64 X 16 ('Table A3 refers) . 





Figure A2.- (Program of Table A2) - Lower half (O^iLz^'i-l = 3) of the extended PHI array showing 
row and column designations of "chunks." IROW 1 and 2 of JCOLUMN 1 and '2 constitute the active 
PHI array. The numbers on "chunks" of JCOLUMN 1 and 2 indicate the numbers of the disk files 
on which those chunks are stored. The "chunks" of JCOLUMN 3 and 4 do not require disk file 
storage. 


o^ 



Storage on disk 
files 1,2,5&6 


"do loop" transfer 


Figure A3.- 
overlays. 



(■program of Table A2) - Arrays dimensioned in the initializing and star advancing 
The numbers on the "chunks" indicate the disk files on which they are stored. 




Figure A4.- (Program of Table A2) - Arrays dimensioned in GETH overlay, which performs a Fourier 
transform of the Greens function H and stores the resulting II- . (Letters A, B, C 

and D are referenced by Table A4.) 













Read from 



(in common) 


Storage on 
disk files 1-8 



RH01 array RH02 array RH03 array RH04 array 


Figure AS.- (Program o£ Table A2) - Arrays dimensioned in the GETPHI overlay, which solves 

foh the potential of an isolated galaxy. 


cn 

KD 




-HI M (f| 


Figure A6.- (Program of Table A2) - Alignment of arrays RHOl, RH02, RH03, and RH04 during calls by 
the GETPHI overlay to its subroutines. Although the active PHI array and the extended PHI array 
are not dimensioned within the GETPHI overlay, their projections on the planes x == 0, y « 0, 
and- z a 0 are represented by dashed and solid lines, respectively. Axes labels represent 
a\ subscripts of array elements which are untransformed C^,y,z), transformed either, 

as appropriate. 



